module fiona_input_mod

  use fiona_common_mod
  use fiona_buffer_mod
  use fiona_dataset_mod
  use fiona_objects_mod

  implicit none

  interface fiona_get_att
    module procedure fiona_get_att_s
    module procedure fiona_get_att_i2
    module procedure fiona_get_att_i4
    module procedure fiona_get_att_i8
    module procedure fiona_get_att_r4
    module procedure fiona_get_att_r8
    module procedure fiona_get_var_att_s
  end interface fiona_get_att

  interface fiona_input
    module procedure fiona_input_s
    module procedure fiona_input_0d_i2
    module procedure fiona_input_0d_i4
    module procedure fiona_input_0d_r4
    module procedure fiona_input_0d_r8
    module procedure fiona_input_1d_i2
    module procedure fiona_input_1d_i4
    module procedure fiona_input_1d_r4
    module procedure fiona_input_1d_r8
    module procedure fiona_input_2d_i2
    module procedure fiona_input_2d_i4
    module procedure fiona_input_2d_r4
    module procedure fiona_input_2d_r8
    module procedure fiona_input_3d_i2
    module procedure fiona_input_3d_i4
    module procedure fiona_input_3d_r4
    module procedure fiona_input_3d_r8
    module procedure fiona_input_4d_i2
    module procedure fiona_input_4d_i4
    module procedure fiona_input_4d_r4
    module procedure fiona_input_4d_r8
    module procedure fiona_input_5d_i2
    module procedure fiona_input_5d_i4
    module procedure fiona_input_5d_r4
    module procedure fiona_input_5d_r8
    module procedure fiona_input_6d_i2
    module procedure fiona_input_6d_i4
    module procedure fiona_input_6d_r4
    module procedure fiona_input_6d_r8
    module procedure fiona_input_7d_i2
    module procedure fiona_input_7d_i4
    module procedure fiona_input_7d_r4
    module procedure fiona_input_7d_r8
  end interface fiona_input

  interface fiona_input_range
    module procedure fiona_input_range_1d_r8
    module procedure fiona_input_range_2d_r8
    module procedure fiona_input_range_3d_r8
  end interface fiona_input_range

contains

  subroutine fiona_open_dataset(dataset_name, file_path, mpi_comm, ngroups, async)

    character(*), intent(in) :: dataset_name
    character(*), intent(in) :: file_path
    integer, intent(in), optional :: mpi_comm
    integer, intent(in), optional :: ngroups
    logical, intent(in), optional :: async

    type(dataset_type), pointer :: dataset
    type(dim_type), pointer :: dim
    type(var_type), pointer :: var
    integer ierr, ncid, ndim, nvar, i, j
    character(NF90_MAX_NAME), allocatable :: dim_names(:), var_names(:), var_dim_names(:,:)
    integer, allocatable :: dim_sizes(:), dimids(:), var_ndim(:)

    if (datasets%hashed(dataset_name)) then
      call log_error('Already open dataset ' // trim(dataset_name) // '!', __FILE__, __LINE__)
    end if

    allocate(dataset)
    call datasets%insert(dataset_name, dataset)
    deallocate(dataset)
    dataset => get_dataset(dataset_name)

    call dataset%init(dataset_name, mpi_comm=mpi_comm, ngroups=ngroups, async=async)
    dataset%file_path = file_path

    ! Inquire dimensions and variables in group root processes.
    if (dataset%is_group_root) then
      ierr = NF90_OPEN(file_path, NF90_NOWRITE, ncid)
      call handle_error(ierr, 'Failed to open NetCDF file ' // trim(file_path) // '!', __FILE__, __LINE__)
      ierr = NF90_INQUIRE(ncid, nDimensions=ndim, nVariables=nvar)
      call handle_error(ierr, 'Failed to inquire number of dimensions and variables in NetCDF file ' // trim(file_path) // '!', __FILE__, __LINE__)
      ! Create dimensions.
      allocate(dim_names(ndim))
      allocate(dim_sizes(ndim))
      do i = 1, ndim
        ierr = NF90_INQUIRE_DIMENSION(ncid, i, name=dim_names(i), len=dim_sizes(i))
        call handle_error(ierr, 'Failed to inquire dimension ' // trim(dim_names(i)) // &
                          ' in NetCDF file ' // trim(file_path) // '!', __FILE__, __LINE__)
        allocate(dim)
        call dataset%dims%insert(dim_names(i), dim)
        deallocate(dim)
        dim => dataset%get_dim(dim_names(i))
        call dim%init(name=dim_names(i), id=i, size=dim_sizes(i))
      end do
      ! Create variables.
      allocate(var_names(nvar))
      allocate(var_ndim(nvar))
      do i = 1, nvar
        ierr = NF90_INQUIRE_VARIABLE(ncid, i, name=var_names(i), ndims=var_ndim(i))
        call handle_error(ierr, 'Failed to inquire variable ' // trim(var_names(i)) // &
                          ' in NetCDF file ' // trim(file_path) // '!', __FILE__, __LINE__)
        allocate(dimids(var_ndim(i)))
        allocate(var_dim_names(var_ndim(i),1))
        ierr = NF90_INQUIRE_VARIABLE(ncid, i, dimids=dimids)
        call handle_error(ierr, 'Failed to inquire variable ' // trim(var_names(i)) // &
                          ' in NetCDF file ' // trim(file_path) // '!', __FILE__, __LINE__)
        do j = 1, var_ndim(i)
          ierr = NF90_INQUIRE_DIMENSION(ncid, dimids(j), name=var_dim_names(j,1))
          call handle_error(ierr, 'Failed to inquire dimension ' // trim(var_dim_names(j,1)) // &
                            ' in NetCDF file ' // trim(file_path) // '!', __FILE__, __LINE__)
        end do
        allocate(var)
        call dataset%vars%insert(var_names(i), var)
        deallocate(var)
        var => dataset%get_var(var_names(i))
        call var%init(name=var_names(i), id=i, ndim=var_ndim(i))
        do j = 1, var_ndim(i)
          var%dims(j)%ptr => dataset%get_dim(var_dim_names(j,1))
        end do
        deallocate(dimids)
        deallocate(var_dim_names)
      end do
      ierr = NF90_CLOSE(ncid)
      ! Record dimension names of each variable for later use.
      allocate(var_dim_names(maxval(var_ndim),nvar))
      do i = 1, nvar
        var => dataset%get_var(var_names(i))
        do j = 1, size(var%dims)
          var_dim_names(j,i) = var%dims(j)%ptr%name
        end do
      end do
    end if

    ! Broadcast dimensions and variables to other processes in the group.
#ifdef HAS_MPI
    if (dataset%group_mpi_comm /= MPI_COMM_NULL) then
      call MPI_BCAST(ndim, 1, MPI_INTEGER, 0, dataset%group_mpi_comm, ierr)
      call MPI_BCAST(nvar, 1, MPI_INTEGER, 0, dataset%group_mpi_comm, ierr)
      if (.not. dataset%is_group_root) then
        allocate(dim_names(ndim))
        allocate(dim_sizes(ndim))
        allocate(var_names(nvar))
        allocate(var_ndim (nvar))
      end if
      call MPI_BCAST(var_ndim, nvar, MPI_INTEGER, 0, dataset%group_mpi_comm, ierr)
      if (.not. dataset%is_group_root) allocate(var_dim_names(maxval(var_ndim),nvar))
      call MPI_BCAST(dim_names, NF90_MAX_NAME * ndim, MPI_CHARACTER, 0, dataset%group_mpi_comm, ierr)
      call MPI_BCAST(dim_sizes, ndim, MPI_INTEGER, 0, dataset%group_mpi_comm, ierr)
      call MPI_BCAST(var_names, NF90_MAX_NAME * nvar, MPI_CHARACTER, 0, dataset%group_mpi_comm, ierr)
      call MPI_BCAST(var_dim_names, NF90_MAX_NAME * maxval(var_ndim) * nvar, MPI_CHARACTER, 0, dataset%group_mpi_comm, ierr)
    end if
#endif

    ! Create dims and vars in the dataset in other processes.
    if (.not. dataset%is_group_root) then
      do i = 1, ndim
        allocate(dim)
        call dataset%dims%insert(dim_names(i), dim)
        deallocate(dim)
        dim => dataset%get_dim(dim_names(i))
        call dim%init(name=dim_names(i), size=dim_sizes(i))
      end do
      do i = 1, nvar
        allocate(var)
        call dataset%vars%insert(var_names(i), var)
        deallocate(var)
        var => dataset%get_var(var_names(i))
        call var%init(name=var_names(i), ndim=var_ndim(i))
        do j = 1, var_ndim(i)
          var%dims(j)%ptr => dataset%get_dim(var_dim_names(j,i))
        end do
      end do
    end if

    deallocate(dim_names)
    deallocate(dim_sizes)
    deallocate(var_names)
    deallocate(var_ndim)
    deallocate(var_dim_names)

  end subroutine fiona_open_dataset

  subroutine fiona_start_input(dataset_name)

    character(*), intent(in) :: dataset_name

    type(dataset_type), pointer :: dataset

    dataset => get_dataset(dataset_name)

    call dataset%open()

  end subroutine fiona_start_input

  subroutine fiona_get_dim(dataset_name, name, size)

    character(*), intent(in) :: dataset_name
    character(*), intent(in) :: name
    integer, intent(out), optional :: size

    type(dataset_type), pointer :: dataset
    type(dim_type), pointer :: dim
    integer ierr, dimid

    dataset => get_dataset(dataset_name)

    if (dataset%dims%hashed(name)) then
      dim => dataset%get_dim(name)
    else
      call log_error('Dimension ' // trim(name) // ' not found in dataset ' // trim(dataset_name) // '!', __FILE__, __LINE__)
    end if

    if (present(size)) size = dim%size

  end subroutine fiona_get_dim

  subroutine fiona_get_att_s(dataset_name, att_name, value)

    character(*), intent(in) :: dataset_name
    character(*), intent(in) :: att_name
    character(*), intent(out) :: value

    type(dataset_type), pointer :: dataset
    integer ierr

    dataset => get_dataset(dataset_name)

    if (dataset%is_group_root) then
      ierr = NF90_GET_ATT(dataset%id, NF90_GLOBAL, att_name, value)
      call handle_error(ierr, 'Failed to get attribute "' // trim(att_name) // &
                        '" from file ' // trim(dataset%file_path) // '!', __FILE__, __LINE__)
    end if
#ifdef HAS_MPI
    if (dataset%group_mpi_comm /= MPI_COMM_NULL) then
      call MPI_BCAST(value, len(value), MPI_CHAR, 0, dataset%group_mpi_comm, ierr)
    end if
#endif

  end subroutine fiona_get_att_s

  subroutine fiona_get_var_att_s(dataset_name, var_name, att_name, value)

    character(*), intent(in ) :: dataset_name
    character(*), intent(in ) :: var_name
    character(*), intent(in ) :: att_name
    character(*), intent(out) :: value

    type(dataset_type), pointer :: dataset
    type(var_type), pointer :: var
    integer ierr

    character(30) name
    integer natt, i

    dataset => get_dataset(dataset_name)

    if (dataset%is_group_root) then
      if (dataset%vars%hashed(var_name)) then
        var => dataset%get_var(var_name)
      else
        allocate(var)
        var%name = var_name
        call dataset%vars%insert(var_name, var)
        deallocate(var)
        var => dataset%get_var(var_name)

        ierr = NF90_INQ_VARID(dataset%id, var_name, var%id)
        call handle_error(ierr, 'Failed to inquire id of variable "' // trim(var_name) // '"!')
      end if

      ierr = NF90_GET_ATT(dataset%id, var%id, att_name, value)
      if (ierr /= NF90_NOERR) then
        ierr = NF90_INQUIRE_VARIABLE(dataset%id, var%id, natts=natt)
        call handle_error(ierr, 'Failed to inquire attribute number of variable "' // &
                          trim(var_name) // '"!', __FILE__, __LINE__)
        do i = 1, natt
          ierr = NF90_INQ_ATTNAME(dataset%id, var%id, i, name)
          call handle_error(ierr, 'Failed to inquire attribute name of variable "' // &
                            trim(var_name) // '"!', __FILE__, __LINE__)
        end do
      end if
      call handle_error(ierr, 'Failed to get attribute "' // trim(att_name) // '" of variable "' // &
                        trim(var_name) // '" from file ' // trim(dataset%file_path) // '!', __FILE__, __LINE__)
    end if
#ifdef HAS_MPI
    if (dataset%group_mpi_comm /= MPI_COMM_NULL) then
      call MPI_BCAST(value, len(value), MPI_CHAR, 0, dataset%group_mpi_comm, ierr)
    end if
#endif

  end subroutine fiona_get_var_att_s

  subroutine fiona_get_att_i2(dataset_name, att_name, value)

    character(*), intent(in) :: dataset_name
    character(*), intent(in) :: att_name
    integer(2), intent(out) :: value

    type(dataset_type), pointer :: dataset
    integer ierr

    dataset => get_dataset(dataset_name)

    if (dataset%is_group_root) then
      ierr = NF90_GET_ATT(dataset%id, NF90_GLOBAL, att_name, value)
      call handle_error(ierr, 'Failed to get global attribute "' // trim(att_name) // &
                        '" from file ' // trim(dataset%file_path) // '!', __FILE__, __LINE__)
    end if
#ifdef HAS_MPI
    if (dataset%group_mpi_comm /= MPI_COMM_NULL) then
      call MPI_BCAST(value, 1, MPI_SHORT, 0, dataset%group_mpi_comm, ierr)
    end if
#endif

 end subroutine fiona_get_att_i2

  subroutine fiona_get_att_i4(dataset_name, att_name, value)

    character(*), intent(in) :: dataset_name
    character(*), intent(in) :: att_name
    integer(4), intent(out) :: value

    type(dataset_type), pointer :: dataset
    integer ierr

    dataset => get_dataset(dataset_name)

    if (dataset%is_group_root) then
      ierr = NF90_GET_ATT(dataset%id, NF90_GLOBAL, att_name, value)
      call handle_error(ierr, 'Failed to get global attribute "' // trim(att_name) // &
                        '" from file ' // trim(dataset%file_path) // '!', __FILE__, __LINE__)
    end if
#ifdef HAS_MPI
    if (dataset%group_mpi_comm /= MPI_COMM_NULL) then
      call MPI_BCAST(value, 1, MPI_INT, 0, dataset%group_mpi_comm, ierr)
    end if
#endif

  end subroutine fiona_get_att_i4

  subroutine fiona_get_att_i8(dataset_name, att_name, value)

    character(*), intent(in) :: dataset_name
    character(*), intent(in) :: att_name
    integer(8), intent(out) :: value

    type(dataset_type), pointer :: dataset
    integer ierr

    dataset => get_dataset(dataset_name)

    if (dataset%is_group_root) then
      ierr = NF90_GET_ATT(dataset%id, NF90_GLOBAL, att_name, value)
      call handle_error(ierr, 'Failed to get global attribute "' // trim(att_name) // &
                        '" from file ' // trim(dataset%file_path) // '!', __FILE__, __LINE__)
    end if
#ifdef HAS_MPI
    if (dataset%group_mpi_comm /= MPI_COMM_NULL) then
      call MPI_BCAST(value, 1, MPI_LONG, 0, dataset%group_mpi_comm, ierr)
    end if
#endif

  end subroutine fiona_get_att_i8

  subroutine fiona_get_att_r4(dataset_name, att_name, value)

    character(*), intent(in) :: dataset_name
    character(*), intent(in) :: att_name
    real(4), intent(out) :: value

    type(dataset_type), pointer :: dataset
    integer ierr

    dataset => get_dataset(dataset_name)

    if (dataset%is_group_root) then
      ierr = NF90_GET_ATT(dataset%id, NF90_GLOBAL, att_name, value)
      call handle_error(ierr, 'Failed to get global attribute "' // trim(att_name) // &
                        '" from file ' // trim(dataset%file_path) // '!', __FILE__, __LINE__)
    end if
#ifdef HAS_MPI
    if (dataset%group_mpi_comm /= MPI_COMM_NULL) then
      call MPI_BCAST(value, 1, MPI_FLOAT, 0, dataset%group_mpi_comm, ierr)
    end if
#endif

  end subroutine fiona_get_att_r4

  subroutine fiona_get_att_r8(dataset_name, att_name, value)

    character(*), intent(in) :: dataset_name
    character(*), intent(in) :: att_name
    real(8), intent(out) :: value

    type(dataset_type), pointer :: dataset
    integer ierr

    dataset => get_dataset(dataset_name)

    if (dataset%is_group_root) then
      ierr = NF90_GET_ATT(dataset%id, NF90_GLOBAL, att_name, value)
      call handle_error(ierr, 'Failed to get global attribute "' // trim(att_name) // &
                        '" from file ' // trim(dataset%file_path) // '!', __FILE__, __LINE__)
    end if
#ifdef HAS_MPI
    if (dataset%group_mpi_comm /= MPI_COMM_NULL) then
      call MPI_BCAST(value, 1, MPI_DOUBLE, 0, dataset%group_mpi_comm, ierr)
    end if
#endif

  end subroutine fiona_get_att_r8

  subroutine fiona_input_s(dataset_name, var_name, value, start, count, fatal)

    character(*), intent(in) :: dataset_name
    character(*), intent(in) :: var_name
    character(*), intent(out) :: value
    integer, intent(in) :: start(2)
    integer, intent(in) :: count(2)
    logical, intent(in), optional :: fatal

    type(dataset_type), pointer :: dataset
    type(var_type), pointer :: var
    logical fatal_opt
    integer ierr, varid

    ! All count should be greater than zero!
    if (any(count <= 0)) then
      call log_error("Variable " // trim(var_name) // " dimension is problematic!")
    end if

    fatal_opt = .false.; if (present(fatal)) fatal_opt = fatal

    dataset => get_dataset(dataset_name)
    var => dataset%get_var(var_name)
    if (dataset%is_group_root) then
      ierr = NF90_INQ_VARID(dataset%id, var_name, varid)
      call handle_error(ierr, 'No variable "' // trim(var_name) // '" in dataset "' // &
                        trim(dataset%file_path) // '"!', __FILE__, __LINE__)
      ierr = NF90_GET_VAR(dataset%id, varid, value, start=start, count=count)
      if (ierr /= NF90_NOERR) then
        if (fatal_opt) then
          call log_error('Failed to read variable "' // trim(var_name) // '" in dataset "' // &
                        trim(dataset%file_path) // '"!', __FILE__, __LINE__)
        else
          call log_warning('Failed to read variable "' // trim(var_name) // '" in dataset "' // &
                        trim(dataset%file_path) // '"!')
        end if
      end if
#ifdef HAS_MPI
      if (dataset%is_parallel) then
        call MPI_BCAST(value, len(value), MPI_CHAR, 0, dataset%group_mpi_comm, ierr)
        if (ierr /= MPI_SUCCESS) then
          if (fatal_opt) then
            call log_error('Failed to broadcast variable "' // trim(var_name) // '" in dataset "' // &
                          trim(dataset%file_path) // '"!', __FILE__, __LINE__)
          else
            call log_warning('Failed to broadcast variable "' // trim(var_name) // '" in dataset "' // &
                          trim(dataset%file_path) // '"!')
          end if
        end if
      end if
#endif
    end if

  end subroutine fiona_input_s

  subroutine fiona_end_input(dataset_name)

    character(*), intent(in) :: dataset_name

    type(dataset_type), pointer :: dataset

    dataset => get_dataset(dataset_name)

    call dataset%close()

    call datasets%remove(trim(dataset_name))

  end subroutine fiona_end_input

#include "fiona_input.F90"

  subroutine read_coord(dataset, dim, x)

    type(dataset_type), intent(in) :: dataset
    type(dim_type), intent(in) :: dim
    real(8), intent(inout), allocatable :: x(:)

    if (allocated(x)) deallocate(x)
    allocate(x(dim%size))
    call fiona_input(dataset%name, dim%name, x)
    if (dim%flip) x = x(dim%size:1:-1)
    if (abs(x(1      ) - dim%span(1)) < 1.0e-12) x(1      ) = dim%span(1)
    if (abs(x(size(x)) - dim%span(2)) < 1.0e-12) x(size(x)) = dim%span(2)

  end subroutine read_coord

  recursive subroutine find_range_indices(dataset, dim, x_in, nx, is, ie, offset, coord_range)

    type(dataset_type), intent(in) :: dataset
    type(dim_type), intent(in) :: dim
    real(8), intent(in) :: x_in(dim%size)
    integer, intent(out) :: nx
    integer, intent(out), allocatable :: is(:)
    integer, intent(out), allocatable :: ie(:)
    real(8), intent(out), allocatable, optional :: offset(:)
    real(8), intent(in), optional :: coord_range(2)

    real(8) x(dim%size)
    real(8) new_coord_range(2)
    real(8) r1, r2, x1, x2
    integer i, p, mode

    x = x_in

    if (present(coord_range)) then
      mode = -1
      x1 = x_in(1); x2 = x_in(size(x_in)); r1 = coord_range(1); r2 = coord_range(2)
      if (x1 <= r1 .and. r2 <= x2) then
          ! o---x-----------x---o
          ! x1  r1          r2  x2
        mode = 1
      else if (r1 <= x1 .and. x1 < r2 .and. r2 <= x2) then
        ! x---o-----------x---o
        ! r1  x1          r2  x2
        mode = 2
      else if (x1 <= r1 .and. r1 < x2 .and. x2 <= r2) then
        ! o---x-----------o---x
        ! x1  r1          x2  r2
        mode = 3
      else if (r1 < x1 .and. x2 < r2) then
          ! x---o-----------o---x
          ! r1  x1          x2  r2
        mode = 4
      else if (dim%cyclic) then
        if (x1 <= r1 - dim%length .and. r1 - dim%length < x2 .and. x2 <= r2 - dim%length) then
          ! o---x-----------o---x
          ! x1  r1          x2  r2
          mode = 5
          x = x + dim%length
        else if (x1 <= r1 - dim%length .and. r1 - dim%length < x2 .and. r2 - dim%length <= x2) then
          ! o---x-----------x---o
          ! x1  r1          r2  x2
          mode = 6
          x = x + dim%length
        end if
      else
        mode = 0
      end if
    else
      mode = 0
    end if
    select case (mode)
    case (0)
      allocate(is(1), ie(1)); is = 0; ie = 0
      if (present(offset)) allocate(offset(1))
      is(1) = 1
      ie(1) = dim%size
      if (present(offset)) offset(1) = 0
    case (1, 6)
      ! o---x-----------x---o
      ! x1  r1          r2  x2
      allocate(is(1), ie(1)); is = 0; ie = 0
      if (present(offset)) allocate(offset(1))
      do i = 1, dim%size
        if (x(i) >= r1) then
          is(1) = max(i - 1, 1)
          exit
        end if
      end do
      do i = is(1) + 1, dim%size
        if (x(i) >= r2) then
          ie(1) = i
          exit
        end if
      end do
      if (ie(1) == 0) ie(1) = dim%size
      if (present(offset)) offset(1) = merge(0.0d0, dim%length, mode == 1)
    case (2)
      if (dim%cyclic) then
        allocate(is(2), ie(2))
      else
        allocate(is(1), ie(1))
      end if
      is = 0; ie = 0
      if (present(offset)) allocate(offset(merge(2, 1, dim%cyclic)))
      if (dim%cyclic) then
        do i = dim%size, 1, -1
          if (x(i) < r1 + dim%length) then
            is(1) = i
            exit
          end if
        end do
        ie(1) = dim%size
        if (present(offset)) offset(1) = -dim%length
        is(2) = 1
        do i = 1, dim%size
          if (x(i) > r2) then
            ie(2) = i
            exit
          end if
        end do
        if (ie(2) == 0) ie(2) = dim%size
        if (present(offset)) offset(2) = 0
      else
        is(1) = 1
        do i = 1, dim%size
          if (x(i) > r2) then
            ie(1) = i
            exit
          end if
        end do
        if (ie(1) == 0) ie(1) = dim%size
        if (present(offset)) offset(1) = 0
      end if
    case (3, 5)
      ! o---x-----------o---x
      ! x1  r1          x2  r2
      if (dim%cyclic) then
        allocate(is(2), ie(2))
      else
        allocate(is(1), ie(1))
      end if
      is = 0; ie = 0
      if (present(offset)) allocate(offset(merge(2, 1, dim%cyclic)))
      if (dim%cyclic) then
        do i = 1, dim%size
          if (x(i) >= r1) then
            is(1) = max(i - 1, 1)
            exit
          end if
        end do
        ie(1) = dim%size
        if (present(offset)) offset(1) = 0
        is(2) = 1
        do i = 1, dim%size
          if (x(i) >= r2 - dim%length) then
            ie(2) = i
            exit
          end if
        end do
        if (ie(2) == 0) ie(2) = 1
        if (present(offset)) offset(2) = dim%length
      else
        do i = 1, dim%size
          if (x(i) > r1) then
            is(1) = i
            exit
          end if
        end do
        ie(1) = dim%size
        if (present(offset)) offset(1) = 0
      end if
    case (4)
      if (dim%cyclic) then
        allocate(is(3), ie(3))
      else
        allocate(is(1), ie(1))
      end if
      is = 0; ie = 0
      if (present(offset)) allocate(offset(merge(3, 1, dim%cyclic)))
      if (dim%cyclic) then
        do i = dim%size, 1, -1
          if (x(i) < r1 + dim%length) then
            is(1) = i
            exit
          end if
        end do
        if (is(1) == 0) is(1) = 1
        ie(1) = dim%size
        if (present(offset)) offset(1) = -dim%length
        is(2) = 1
        ie(2) = dim%size
        if (present(offset)) offset(2) = 0
        is(3) = 1
        do i = 1, dim%size
          if (x(i) > r2 - dim%length) then
            ie(3) = i
            exit
          end if
        end do
        if (ie(3) == 0) ie(3) = dim%size
        if (present(offset)) offset(3) = dim%length
      else
        is(1) = 1
        ie(1) = dim%size
        if (present(offset)) offset(1) = 0
      end if
    end select
    if (dataset%async .and. dataset%is_group_root) then
      if (allocated(is)) deallocate(is, ie)
      allocate(is(1)); is = 1
      allocate(ie(1)); ie = 1
      nx = 1
    else
      if (dim%flip) then
        is = dim%size - is + 1
        ie = dim%size - ie + 1
        do p = 1, size(is)
          i = is(p)
          is(p) = ie(p)
          ie(p) = i
        end do
      end if
      nx = sum(ie - is) + size(is)
    end if

    if (.not. all(is < ie)) then
      if (present(coord_range)) then
        new_coord_range(1) = coord_range(1) - 0.1d0 * (coord_range(2) - coord_range(1))
        new_coord_range(2) = coord_range(2) + 0.1d0 * (coord_range(2) - coord_range(1))
        deallocate(is, ie)
        if (present(offset)) deallocate(offset)
        call find_range_indices(dataset, dim, x_in, nx, is, ie, offset=offset, coord_range=new_coord_range)
      else
        write(*, *) '[Error]: Unexpected error in find_range_indices!'
        stop 1
      end if
    end if

  end subroutine find_range_indices

  subroutine fiona_input_range_1d_r8(dataset_name, var_name, array, coord_range)

    character(*), intent(in) :: dataset_name
    character(*), intent(in) :: var_name
    real(8), intent(inout), allocatable :: array(:)
    real(8), intent(in) :: coord_range(2)

    type(dataset_type), pointer :: dataset
    type(var_type), pointer :: var
    type(dim_type), pointer :: dim
    real(8), allocatable :: x(:)
    real(8), allocatable :: offset(:)
    integer, allocatable :: is(:), ie(:)
    integer i, nx, i1, i2, max_is_size, ierr
    real(8), allocatable :: tmp(:)

    dataset => get_dataset(dataset_name)
    var => dataset%get_var(var_name)
    dim => var%dims(1)%ptr

    call read_coord(dataset, dim, x)
    if (var_name == dim%name) then
      call find_range_indices(dataset, dim, x, nx, is, ie, offset=offset, coord_range=coord_range)
    else
      call find_range_indices(dataset, dim, x, nx, is, ie, coord_range=coord_range)
    end if
    deallocate(x)
    if (allocated(array)) deallocate(array)
    allocate(array(nx))

#ifdef HAS_MPI
    call MPI_ALLREDUCE(size(is), max_is_size, 1, MPI_INTEGER, MPI_MAX, dataset%global_mpi_comm, ierr)
#endif

    allocate(tmp(1))
    i1 = 1
    do i = 1, size(is)
      i2 = i1 + ie(i) - is(i)
      call fiona_input(dataset_name, var_name, array(i1:i2), start=[is(i)], count=[ie(i)-is(i)+1])
      if (var_name == dim%name .and. .not. (dataset%async .and. dataset%is_group_root)) then
        if (offset(i) /= 0) array(i1:i2) = array(i1:i2) + offset(i)
      end if
      i1 = i2 + 1
    end do
    do i = size(is) + 1, max_is_size
      call fiona_input(dataset_name, var_name, tmp, start=[1], count=[1])
    end do
    deallocate(tmp)
    if (dim%flip) then
      allocate(tmp(nx))
      tmp = array(nx:1:-1)
      array = tmp
      deallocate(tmp)
    end if

    if (allocated(is)) deallocate(is, ie)
    if (allocated(offset)) deallocate(offset)

  end subroutine fiona_input_range_1d_r8

  subroutine fiona_input_range_2d_r8(dataset_name, var_name, array, coord_range_1, coord_range_2)

    character(*), intent(in) :: dataset_name
    character(*), intent(in) :: var_name
    real(8), intent(inout), allocatable :: array(:,:)
    real(8), intent(in), optional :: coord_range_1(2)
    real(8), intent(in), optional :: coord_range_2(2)

    type(dataset_type), pointer :: dataset
    type(var_type), pointer :: var
    type(dim_type), pointer :: dim1, dim2
    real(8), allocatable :: x1(:), x2(:)
    integer, allocatable :: is(:), ie(:), js(:), je(:)
    integer i, j, nx1, nx2, i1, i2, j1, j2, max_is_size, max_js_size, ierr
    real(8), allocatable :: tmp(:,:)

    dataset => get_dataset(dataset_name)
    var => dataset%get_var(var_name)
    dim1 => var%dims(1)%ptr
    dim2 => var%dims(2)%ptr

    call read_coord(dataset, dim1, x1)
    call find_range_indices(dataset, dim1, x1, nx1, is, ie, coord_range=coord_range_1)
    call read_coord(dataset, dim2, x2)
    call find_range_indices(dataset, dim2, x2, nx2, js, je, coord_range=coord_range_2)
    deallocate(x1, x2)
    if (allocated(array)) deallocate(array)
    allocate(array(nx1,nx2))

#ifdef HAS_MPI
    call MPI_ALLREDUCE(size(is), max_is_size, 1, MPI_INTEGER, MPI_MAX, dataset%global_mpi_comm, ierr)
    call MPI_ALLREDUCE(size(js), max_js_size, 1, MPI_INTEGER, MPI_MAX, dataset%global_mpi_comm, ierr)
#endif

    allocate(tmp(1,1))
    i1 = 1
    do i = 1, size(is)
      i2 = i1 + ie(i) - is(i)
      j1 = 1
      do j = 1, size(js)
        j2 = j1 + je(j) - js(j)
        call fiona_input(dataset_name, var_name, array(i1:i2,j1:j2), start=[is(i),js(j)], count=[ie(i)-is(i)+1,je(j)-js(j)+1])
        j1 = j2 + 1
      end do
      do j = size(js) + 1, max_js_size
        call fiona_input(dataset_name, var_name, tmp, start=[1,1], count=[1,1])
      end do
      i1 = i2 + 1
    end do
    do i = size(is) + 1, max_is_size
      do j = 1, max_js_size
        call fiona_input(dataset_name, var_name, tmp, start=[1,1], count=[1,1])
      end do
    end do
    deallocate(tmp)
    if (dim1%flip .or. dim2%flip) allocate(tmp(nx1,nx2))
    if (dim1%flip) then
      tmp = array(nx1:1:-1,:)
      array = tmp
    end if
    if (dim2%flip) then
      tmp = array(:,nx2:1:-1)
      array = tmp
    end if
    if (allocated(tmp)) deallocate(tmp)

    if (allocated(is)) deallocate(is, ie, js, je)

  end subroutine fiona_input_range_2d_r8

  subroutine fiona_input_range_3d_r8(dataset_name, var_name, array, coord_range_1, coord_range_2, coord_range_3)

    character(*), intent(in) :: dataset_name
    character(*), intent(in) :: var_name
    real(8), intent(inout), allocatable :: array(:,:,:)
    real(8), intent(in), optional :: coord_range_1(2)
    real(8), intent(in), optional :: coord_range_2(2)
    real(8), intent(in), optional :: coord_range_3(2)

    type(dataset_type), pointer :: dataset
    type(var_type), pointer :: var
    type(dim_type), pointer :: dim1, dim2, dim3
    real(8), allocatable :: x1(:), x2(:), x3(:)
    integer, allocatable :: is(:), ie(:), js(:), je(:), ks(:), ke(:)
    integer i, j, k, nx1, nx2, nx3, i1, i2, j1, j2, k1, k2, max_is_size, max_js_size, max_ks_size, ierr
    real(8), allocatable :: tmp(:,:,:)

    dataset => get_dataset(dataset_name)
    var => dataset%get_var(var_name)
    dim1 => var%dims(1)%ptr
    dim2 => var%dims(2)%ptr
    dim3 => var%dims(3)%ptr

    call read_coord(dataset, dim1, x1)
    call find_range_indices(dataset, dim1, x1, nx1, is, ie, coord_range=coord_range_1)
    call read_coord(dataset, dim2, x2)
    call find_range_indices(dataset, dim2, x2, nx2, js, je, coord_range=coord_range_2)
    call read_coord(dataset, dim3, x3)
    call find_range_indices(dataset, dim3, x3, nx3, ks, ke, coord_range=coord_range_3)
    deallocate(x1, x2, x3)
    if (allocated(array)) deallocate(array)
    allocate(array(nx1,nx2,nx3))

#ifdef HAS_MPI
    call MPI_ALLREDUCE(size(is), max_is_size, 1, MPI_INTEGER, MPI_MAX, dataset%global_mpi_comm, ierr)
    call MPI_ALLREDUCE(size(js), max_js_size, 1, MPI_INTEGER, MPI_MAX, dataset%global_mpi_comm, ierr)
    call MPI_ALLREDUCE(size(ks), max_ks_size, 1, MPI_INTEGER, MPI_MAX, dataset%global_mpi_comm, ierr)
#endif

    allocate(tmp(1,1,1))
    i1 = 1
    do i = 1, size(is)
      i2 = i1 + ie(i) - is(i)
      j1 = 1
      do j = 1, size(js)
        j2 = j1 + je(j) - js(j)
        k1 = 1
        do k = 1, size(ks)
          k2 = k1 + ke(k) - ks(k)
          call fiona_input(dataset_name, var_name, array(i1:i2,j1:j2,k1:k2), start=[is(i),js(j),ks(k)], count=[ie(i)-is(i)+1,je(j)-js(j)+1,ke(k)-ks(k)+1])
          k1 = k2 + 1
        end do
        do k = size(ks) + 1, max_ks_size
          call fiona_input(dataset_name, var_name, tmp, start=[1,1,1], count=[1,1,1])
        end do
        j1 = j2 + 1
      end do
      do j = size(js) + 1, max_js_size
        do k = 1, max_ks_size
          call fiona_input(dataset_name, var_name, tmp, start=[1,1,1], count=[1,1,1])
        end do
      end do
      i1 = i2 + 1
    end do
    do i = size(is) + 1, max_is_size
      do j = 1, max_js_size
        do k = 1, max_ks_size
          call fiona_input(dataset_name, var_name, tmp, start=[1,1,1], count=[1,1,1])
        end do
      end do
    end do
    deallocate(tmp)
    if (dim1%flip .or. dim2%flip .or. dim3%flip) allocate(tmp(nx1,nx2,nx3))
    if (dim1%flip) then
      tmp = array(nx1:1:-1,:,:)
      array = tmp
    end if
    if (dim2%flip) then
      tmp = array(:,nx2:1:-1,:)
      array = tmp
    end if
    if (dim3%flip) then
      tmp = array(:,:,nx3:1:-1)
      array = tmp
    end if
    if (allocated(tmp)) deallocate(tmp)

    if (allocated(is)) deallocate(is, ie, js, je, ks, ke)

  end subroutine fiona_input_range_3d_r8

end module fiona_input_mod
